function [out] = my_resampleMatrix(Mat,res1,res2,opt) 
	% out = my_resampleMatrix(Mat,res1,res2)
	% -------- custom function ------
	% description:
	%	resample matrix acoording to bigger resolution,for example: from 1 to 
	%	0.05 degree
	% input:
	%	Mat, old matrix, support 3d matrix that 3rd dimension is time
	%		or lat,lon of old matrix
	%	res1, old resolution
	%	res2, new resolution
	%	opt, 'lat' or 'lon' or 'sub' or 'ind', convert index or coordinates
	% output:
	%	out, matrix of new resolution
	% Author: Sid Chen
	% Date: 2022/2/22
	% Log:
	% 	2022/2/22	ver1.0	create
	arguments
		Mat
		res1
		res2
		opt.VarName {mustBeMember(opt.VarName,{'lat','lon','sub','ind'})} = 'lat'
	end

	Scale = res1 ./ res2;
	switch ndims(Mat)
		case 2
			if isVector(Mat) && (opt.VarName == "lat" || opt.VarName == "lon")
				% if is lat/lon
				out = resampleCoor(Mat,Scale);
            elseif isVector(Mat) && (opt.VarName == "sub" || opt.VarName == "ind")
                % if is linear or sub index
                out = getIndex(Mat,Scale);
			else
				% if 2d matrix
				out = resample2D(Mat,Scale);
			end
		case 3
			% if 3rd dimension is time
			Dim3rd = size(Mat,3);
			Size = [size(Mat,[1,2]) .* (Scale), Dim3rd];
			out = zeros(Size);
			for ii = 1:size(Mat,3)
                Sample = squeeze(Mat(:,:,ii));
				out(:,:,ii) = resample2D(Sample,Scale);
			end
		otherwise
			error("dim of Input matrix must be 2 or 3");
	end
end

function out = isVector(Input)
	if any(size(Input) == 1)
		out =  true;
	else
		out = false;
	end
end

function out = resample2D(Mat,Scale)
	Size = size(Mat);
	x = 1:Size(2);
	y = 1:Size(1);
	x_new = getIndex(x,Scale);
	y_new = getIndex(y,Scale);

	[X,Y] = meshgrid(x_new,y_new);
	Index = sub2ind(Size,Y(:),X(:));

	Size_new = Size .* Scale;
	out = reshape(Mat(Index),Size_new);
end

function out = getIndex(index,Scale)
	index = transpose(index(:));
	index2 = repmat(index,Scale,1);
	out = index2(:);
end

function out = resampleCoor(coor,Scale)
	Step = abs(mode(diff(coor)));
	Step_new = Step ./ Scale;
	Sign = sign(coor(end) - coor(1));
	Edge1 = coor(1) - Sign .* Step ./ 2;
	Edge2 = coor(end) + Sign .* Step ./2;
	new_start = Edge1 + Sign .* Step_new ./2;
	new_end = Edge2 - Sign .* Step_new ./2;
	new_num = numel(coor) .* Scale;
	out = linspace(new_start,new_end,new_num);
    out = out(:);
end
